% This script rescales and plots impulse response functions for 
% contemporaneous or news TFP shocks.

% Christoph Gortz


close all
clear all

%% Initialisation

 
shockPlot = 'z';             % Shock process to be plotted
exo_shockPlot = 'ez8';        % Exogenous shock to be plotted
shockOther = 'v';             % Other Shock process with trend
exo_shockOther = 'ev';        % Other Exogenous shock with trend
vars = [                    % Variables to be plotted
    'y      ';
    'c      ';
    'w      ';
    'pi     ';
    'wMarkup';
%     'sc     ';
%     'si     ';
%     'spreadc';
    %
    'inve   ';
    'invec  ';
    'invei  ';
%     'nc     ';
%     'ni     ';
%     'spreadi';
    'qc     ';
    'mcc    ';
    %
    'lab    ';
    'labc   ';
    'labi   ';
%     'lrc    ';
%     'lri    ';
%     'spkc   ';
    %
%     'pinfc  ';
%     'pinfi  ';
%     'w      ';
    'qi     ';
    'mci    ';
    ];
nVars = size(vars,1);
IRFtitle = [                % Title of figures corresponding to variables
    'Output                   ';
    'Consumption              ';
    'Real Wage                ';
%     'C-Sector Financial Claims';
%     'I-Sector Financial Claims';
%     'C-Sector Spread          ';
    'Rel. Price of Investment ';
    'Wage Markup              ';
    %
    'Total Investment         ';
    'C-Sector Investment      ';
    'I-Sector Investment      ';
%     'C-Sector Bank Capital    ';
%     'I-Sector Bank Capital    ';
%     'I-Sector Spread          ';
    'C-Sector Price of Capital';
    'C-Sector Price Markup    ';
    %
    'Total Hours              ';
    'C-Sector Hours           ';
    'I-Sector Hours           ';
%     'C-Sector Leverage Ratio  ';
%     'I-Sector Leverage Ratio  ';
%     'C-Sector Valuation Shock ';
    %
%     'C-Sector Inflation       ';
%     'I-Sector Inflation       ';
%     'Real Wage                ';
%     'C-Sector Price of Capital';
    'I-Sector Price of Capital';
    'I-Sector Price Markup    ';
    ];



% Names of .mat files with IRFs to be plotted in same figure
dataName = ['Output_Baseline1               '       
            'Output_Baseline_noWandPRigidity' % model with no wage and price rigidities
            ]; 
figName = ['IRF_baseline_vs_noWandPRigidity_' , exo_shockPlot];



figDimensions = [3, 5];     % Dimensions of figure - needs to be consistent
                            % with nVars
irfLength = 40;             % Length of generated IRFs 


%% LOAD INFO AND DETREND
for iData = 1:size(dataName,1)
    eval([ 'load  ' dataName(iData,:) ';' ])

% Load variables' irf data into cell
for iVars = 1:nVars
    eval([ 'irfMean(iVars,1) = {oo_.irfs.' strcat(vars(iVars,:)) '_' exo_shockPlot '};' ])
end

for i = 1:numel(irfMean);
    Variables{i,1} = [irfMean{i,1}];
end
% Load shock's irf data
eval([ shockPlot 'Mean = oo_.irfs.' shockPlot '_' exo_shockPlot ';' ])
eval([ shockOther 'Mean = oo_.irfs.' shockOther '_' exo_shockPlot ';' ])



% Transform series of shock z to isolate zd:
vd =[]; vp = 0;
zd =[]; zp = 0;

for m = 1:numel(zMean)
    % Change of v and z in percent
    vp = vMean(1,m)+vp;
    vd = [vd; vp];
    zp = zMean(1,m)+zp;
    zd = [zd; zp];
end    
vd = vd';
zd = zd';


% Descale stationarised variables: K_t = K_t^~*Z_t
yd = Variables{1,1} + (alfac/(1-alfai))*vd + zd;
cd = Variables{2,1} + (alfac/(1-alfai))*vd + zd;
wd = Variables{3,1} + (alfac/(1-alfai))*vd + zd;
pid = Variables{4,1} + zd - ((1-alfac)/(1-alfai))*vd;
wMarkupd = Variables{5,1};
inved = Variables{6,1} + (1/(1-alfai))*vd;
invecd = Variables{7,1} + (1/(1-alfai))*vd;
inveid = Variables{8,1} + (1/(1-alfai))*vd;
qcd = Variables{9,1} - ((1-alfac)/(1-alfai))*vd + zd;
mccd = ( Variables{10,1} )*(-1);  % Inverse as we want to plot the markup instead of marginal cost
labd = Variables{11,1};
labcd = Variables{12,1};
labid = Variables{13,1};
qid = Variables{14,1} - ((1-alfac)/(1-alfai))*vd + zd;
mcid = ( Variables{15,1} )*(-1);  % Inverse as we want to plot the markup instead of marginal cost


% Final Detrended Time Series of IRFs
eval([ 'VariablesDetrended_' dataName(iData,:) '= [yd; cd; wd; pid; wMarkupd; inved; invecd; inveid; qcd; mccd; labd; labcd; labid; qid; mcid];' ]);

end
%% Plot Detrended IRFs

%% Plot IRFs
figure1 = figure('Color',[1 1 1]); % Generates white background around plot

    for iVars = 1:nVars
        varStr = strcat(vars(iVars,:));
        plotVar = strcat(vars(iVars,:));
        irfTitle = strcat(IRFtitle(iVars,:));
        % Subplots
        subplot(figDimensions(1),figDimensions(2),iVars);
        for iData = 1:size(dataName,1) % Loop over different datasets of model versions
        eval([ 'irfData = VariablesDetrended_' dataName(1,:) '(iVars,:);' ]);
        eval([ 'plot(irfData,''-k'',''linewidth'',1.5);' ]);
        hold on
        eval([ 'irfData = VariablesDetrended_' dataName(2,:) '(iVars,:);' ]);
        eval([ 'plot(irfData,''--+b'',''linewidth'',1.5);' ]);
        if iData == 3
            eval([ 'irfData = VariablesDetrended_' dataName(3,:) '(iVars,:);' ]);
            eval([ 'plot(irfData,''-or'',''linewidth'',1.5);' ]);
        end
        end
        axis tight
        hold off
        eval( 'title( irfTitle );' );
    end
set(gcf, 'Position', get(0,'Screensize')); % Maximise Figure
